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ABSTRACT 


Nearshore  depths  for  Waimanalo  Beach,  HI,  are  extracted  from  optical  imagery,  taken  by 
the  WorldView-2  satellite  on  31  March  2011,  by  means  of  automated  wave  kinematics 
bathymetry  (WKB).  Two  sets  of  three  sequential  images  taken  at  intervals  of  about  10 
seconds  are  used  for  the  analyses  herein.  Water  depths  are  calculated  using  a  computer 
program  that  registers  the  images,  estimates  the  currents,  and  then  uses  the  linear 
dispersion  relationship  for  surface  gravity  waves  to  estimate  depth.  Depths  are  generated 
from  close  to  shore  out  to  about  20  meters  depth.  Comparisons  with  SHOALS  LIDAR 
bathymetry  values  show  WKB  depths  are  accurate  to  about  half  a  meter,  with  R2  values 
of  90%,  and  are  frequently  in  the  range  of  10-20  percent  relative  error  for  depths  ranging 
from  2-16  meters. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION . 1 

A.  PURPOSE  OF  RESEARCH . 1 

B.  SPECIFIC  OBJECTIVES . 3 

II.  BACKGROUND . 5 

A.  INTRODUCTION . 5 

B.  THEORY . 5 

C.  SOME  MODERN  METHODS  FOR  DETERMINING  NEARSHORE 

BATHYMETRY . 7 

1.  Airborne  Passive  Optical  System  (Using  Linear  Dispersion 

Theory) . 7 

2.  Synthetic  Aperture  Radar  (Using  Wave  Refraction) . 11 

3.  Aerial  Photography  (Using  Radiometric  Techniques) . 13 

III.  EXPERIMENT  DESIGN . 17 

A.  PROBLEM  DEFINITION . 17 

B.  MATERIALS . 19 

1.  WorldView-2  Satellite . 19 

2.  SHOALS  Bathymetry . 20 

3.  WKB  Algorithm . 20 

C.  METHOD . 21 

1.  Pre-Processing . 24 

2.  WKB  Extraction . 27 

IV.  ANALYSIS  RESULTS . 35 

A.  GENERAL  APPROACH . 35 

B.  IMAGE  SET  1 . 38 

C.  IMAGE  SET  2 . 44 

V.  CONCLUSIONS  AND  RECOMMENDATIONS . 49 

A.  CONCLUSIONS . 49 

B.  RECOMMENDATIONS  FOR  FUTURE  WORK . 50 

APPENDIX  A.  IMAGE  SET  1 . 53 

CASE:  2MSI . 53 

CASE:  3MSI . 58 

CASE:  2PAN . 63 

CASE:  3PAN . 68 

APPENDIX  B.  IMAGE  SET  2 . 73 

CASE:  2MSI . 73 

CASE:  3MSI . 78 

CASE:  2PAN . 83 

CASE:  3PAN . 88 

LIST  OF  REFERENCES . 93 

INITIAL  DISTRIBUTION  LIST . 95 

vii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


LIST  OF  FIGURES 


Figure  1.  Compact  Hydrographic  Airborne  Rapid  Total  Survey  (From  Joint  Airborne 

LIDAR  Bathymetry  Technical  Center  of  Expertise,  2011) . 2 

Figure  2.  Close-up  photo  of  camera  and  turret  as  mounted  on  the  nose  cone  of  the 

Pelican  aircraft  (From  Dugan  et  ah,  2001a) . 8 


Figure  3.  Frequency- wavenumber  slice,  from  Outer  Banks  study,  through  the  power 

spectrum,  oriented  in  the  direction  of  the  primary  swell  waves.  The 
theoretical  dispersion  surface  that  includes  no  current  and  is  for  infinite 
depth  is  represented  by  the  dashed  curve.  The  depth  and  current  are 
determined  by  calculating  a  dispersion  surface  that  fits  the  measured  wave 

spectrum,  denoted  by  the  solid  curve  (From  Dugan  et  al.,  2001b) . 9 

Figure  4.  Frequency- wavenumber  slice,  from  Monterey  Bay  study,  through  the  power 

spectrum,  in  the  direction  of  the  wind.  The  solid  curve  is  the  intersection  of 
the  deep  water  dispersion  surface  (i.e.,  no  current  and  infinite  depth).  The 
wave  energy  is  closely  distributed  along  this  curve,  indicating  that  no 
significant  surface  currents  were  present  during  the  collection  (From  Dugan 


et  al.,  2001a) . 10 

Figure  5.  RMS  error  as  a  function  of  dwell  (From  Dugan  et  al.,  2002) . 1 1 


Figure  6.  Two-scale  model  illustrating  Bragg  waves  (the  smaller  scale  waves) 
embedded  in  and  riding  on  the  large-scale  surface  waves.  Flat  plates,  a  few 
Bragg  wavelengths  in  size,  model  the  radar’s  interaction  with  the  surface. 
The  plates  tilt  and  move  based  on  the  local  slope  and  motion  of  the  large- 


scale  wave  surface  (From  Wackerman  &  Clemente-Colon,  2004) . 12 

Figure  7.  Depth  error  comparison  of  the  ratio  and  modified  methods  (From 

Lyzenga,  1978) . 14 

Figure  8.  WorldView-2  sensor  bands  showing  their  relative  positions  and  overlap  in 

the  electromagnetic  spectrum  (From  DigitalGlobe,  2011b) . 18 

Figure  9.  STK  snap  shot  of  WorldView-2  collection  pass . 21 

Figure  10.  SHOALS  2000  survey  of  Oahu  (From  JALBTCX,  2011) . 22 

Figure  11.  Example  WorldView-2  multispectral  image . 23 

Figure  12.  Example  cropped  area  (in  green) . 24 

Figure  13.  Registration  GUI . 25 

Figure  14.  Deep  water  analysis . 26 

Figure  15.  True  color  image  (left)  with  corresponding  land  mask  (right) . 27 

Figure  16.  Image  Set  1  case  3MSI  WKB  output  showing  wave  direction  (top  left),  a 

true  color  image  (top  right),  extracted  bathymetry  (bottom  left),  and 

extracted  ocean  currents  (bottom  right) . 33 

Figure  17.  Comparison  maps  of  SHOALS  (left)  and  3MSI  WKB  for  Image  Set  1 

(right) . 36 

Figure  18.  Transect  from  Image  Set  1  case  3MSI  showing  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 
and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (bot) . 37 


IX 


Figure  19.  Scatter  plot  from  Image  Set  1  case  3MSI  showing  thinned  data  for  all  depths 

(top)  and  for  just  depths  less  than  15  m  (bottom) . 37 

Figure  20.  Bar  graph  from  Image  Set  1  case  3MSI  showing  the  mean  depth  error  for 

several  depth  bins . 38 

Figure  2 1 .  Mean  depth  errors  for  all  four  Image  Set  1  cases . 39 

Figure  22.  Energy  spectrum  from  Mokapu  Point  buoy  for  March  2011,  showing  larger 
swell  events  earlier  in  the  month,  but  relatively  little  swell  on  the  date  of 

collection  (31  March)  (From  CDIP,  2012) . 42 

Figure  23.  Energy  spectra  from  Mokapu  Point  buoy  for  the  time  of  collection  (left)  and 
a  swell  event  a  few  days  prior  (right),  showing  an  order  of  magnitude 

difference  in  energy  at  the  low  frequencies  (From  CDIP,  2012) . 43 

Figure  24.  Mean  depth  errors  for  all  four  Image  Set  2  cases . 45 

Figure  25.  2MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color  image 
(top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 53 

Figure  26.  Comparison  maps  of  SHOALS  (left)  and  2MSI  WKB  (right) . 54 

Figure  27.  2MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 54 

Figure  28.  2MSI  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 55 

Figure  29.  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 56 

Figure  30.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 56 

Figure  31.  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 57 

Figure  32.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 57 

Figure  33.  3MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color  image 
(top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 58 

Figure  34.  Comparison  maps  of  SHOALS  (left)  and  3MSI  WKB  (right) . 59 

Figure  35.  3MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 59 

Figure  36.  3MSI  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 60 

Figure  37.  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 61 

Figure  38.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 61 

Figure  39.  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 62 


x 


Figure  40.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 62 

Figure  41.  2Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color  image  (top 
right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 63 

Figure  42.  Comparison  maps  of  SHOALS  (left)  and  2Pan  WKB  (right) . 64 

Figure  43.  2Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 64 

Figure  44.  2Pan  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 65 

Figure  45.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 66 

Figure  46.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 66 

Figure  47.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 67 

Figure  48.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 67 

Figure  49.  3Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color  image  (top 
right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 68 

Figure  50.  Comparison  maps  of  SHOALS  (left)  and  3Pan  WKB  (right) . 69 

Figure  5 1 .  3Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 69 

Figure  52.  3Pan  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 70 

Figure  53.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 71 

Figure  54.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 71 

Figure  55.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 72 

Figure  56.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 72 

Figure  57.  2MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color  image 
(top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 73 

Figure  58.  Comparison  maps  of  SHOALS  (left)  and  2MSI  WKB  (right) . 74 

Figure  59.  2MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 74 

Figure  60.  2MSI  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 75 


xi 


Figure  61.  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 76 

Figure  62.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (bot) . 76 

Figure  63.  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 77 

Figure  64.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (bot) . 77 

Figure  65.  3MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color  image 
(top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 78 

Figure  66.  Comparison  maps  of  SHOALS  (left)  and  3MSI  WKB  (right) . 79 

Figure  67.  3MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 79 

Figure  68.  3MSI  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 80 

Figure  69.  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 81 

Figure  70.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (bot) . 81 

Figure  71.  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for  more 

explanation  of  the  panels) . 82 

Figure  72.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (bot) . 82 

Figure  73.  2Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color  image  (top 
right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 83 

Figure  74.  Comparison  maps  of  SHOALS  (left)  and  2Pan  WKB  (right) . 84 

Figure  75.  2Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 84 

Figure  76.  2Pan  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 85 

Figure  77.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 86 

Figure  78.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (bot) . 86 

Figure  79.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 87 

Figure  80.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 
and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (bot) . 87 

xii 


Figure  81.  3Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color  image  (top 
right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean  currents 

(bottom  right) . 88 

Figure  82.  Comparison  maps  of  SHOALS  (left)  and  3Pan  WKB  (right) . 89 

Figure  83.  3Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for  just 

depths  less  than  15  m  (bottom) . 89 

Figure  84.  3Pan  bar  graph  showing  the  mean  depth  error  for  several  depth  bins . 90 

Figure  85.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 91 

Figure  86.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 

and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 91 

Figure  87.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for  more 

explanation  of  the  panels) . 92 

Figure  88.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted  with 
SHOALS  depth  (top);  the  difference  between  them,  or  depth  error  (mid); 
and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) . 92 


xiii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xiv 


LIST  OF  TABLES 


Table  1.  Comparison  of  bathymetry  methods  using  remote  sensing 

(From  Abileah,  2006) . 17 

Table  2.  Key  WorldView-2  specifications  (After  DigitalGlobe,  2011a) . 19 

Table  3.  Summary  of  image  files . 23 

Table  4.  Mean  depth  error  and  relative  error  ranges  for  each  depth  bin  in  Image  Set  1  ...40 

Table  5.  R 2  values  for  Image  Set  1 . 40 

Table  6.  Mean  depth  error  and  relative  error  ranges  for  each  depth  bin  in  Image  Set  2. ..46 
Table  7.  R 2  values  for  Image  Set  2 . 46 


xv 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xvi 


LIST  OF  ACRONYMS  AND  ABBREVIATIONS 


2MSI  Two-image  Multispectral 

2Pan  Two-image  Panchromatic 

3MSI  Three-image  Multispectral 

3Pan  Three-image  Panchromatic 

AROSS  Airborne  Remote  Optical  Spotlight  System 

CDIP  Coastal  Data  Infonnation  Program 

CHARTS  Compact  Hydrographic  Airborne  Rapid  Total  Survey 

GUI  Graphical  User  Interface 

IHO  International  Hydrographic  Organization 

JALBTCX  Joint  Airborne  LIDAR  Bathymetry  Technical  Center  of  Expertise 

LIDAR  Light  Detection  and  Ranging 

MLLW  Mean  Lower  Low  Water 

RMS  Root  Mean  Squared 

SAR  Synthetic  Aperture  Radar 

SHOALS  Scanning  Hydrographic  Operational  Airborne  LIDAR  Survey 

SNR  Signal-to-Noise  Ratio 

STK  Satellite  Tool  Kit 

UTM  Universal  Transverse  Mercator 


WKB  Wave  Kinematics  Bathymetry 


XVII 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


ACKNOWLEDGMENTS 


In  addition  to  my  advisors,  several  people  deserve  many  thanks  for  their  roles  in 
helping  with  the  creation  of  this  thesis.  Mr.  Ron  Abileah  of  jOmegak,  the  inventor  of  the 
algorithm  tested,  provided  a  wealth  of  knowledge  about  the  methods  used  and  was 
extremely  patient  with  the  many  questions  I  had  regarding  the  process  and  the  algorithm. 
His  detailed  explanations  greatly  enhance  my  understanding  and  ability  to  convey  the 
necessary  information  to  the  reader.  Mike  Cook  of  the  NPS  Oceanography  Dept, 
provided  valuable  assistance  with  writing  the  MATLAB  analysis  and  comparison 
program  I  used  to  compare  the  results  of  the  algorithm  to  ground  truth.  Angela  Kim  and 
Krista  Lee  of  the  NPS  Remote  Sensing  Center  were  helpful  in  getting  me  started  and 
pointing  me  in  the  right  direction  when  I  needed  it.  Sarah  Carlisle,  also  of  the  NPS 
Remote  Sensing  Center,  performed  the  tedious  job  of  converting  and  orthorectifying  the 
images,  several  times.  Thank  you  all  for  your  contributions. 

I  would  also  like  to  thank  my  wife,  Katie,  and  my  three  children  for  their  patience 
and  understanding  during  the  four  months  we  were  separated  while  I  devoted  myself  to 
this  project. 


xix 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xx 


I.  INTRODUCTION 


A.  PURPOSE  OF  RESEARCH 

Characterizing  the  environmental  parameters  of  the  battlespace  ahead  of  a 
military  operation  can  enable  the  planners  and  decision  makers  to  understand  the 
potential  impacts  the  environment  might  have.  This  can  help  them  conduct  safer  and 
more  efficient  operations.  One  parameter  that  is  extremely  important  for  certain 
operations  conducted  in  the  nearshore  is  bathymetry.  Amphibious  landings,  mine 
warfare  operations,  reconnaissance  missions,  and  other  special  operations  missions 
performed  in  the  nearshore  region  require  accurate  and  up-to-date  bathymetric 
information.  Without  it,  the  operation  could  be  hampered  by  difficulties,  such  as  vessel 
groundings,  causing  equipment  damage,  personnel  injuries  or  death,  delays,  and  perhaps 
mission  failure. 

Several  methods  and  techniques  are  currently  employed  to  detennine  nearshore 
bathymetry.  Hydrographic  surveys,  especially  those  that  require  International 
Hydrographic  Organization  (IHO)  standards,  such  as  for  navigation  charts,  are  typically 
performed  using  acoustic  systems,  such  as  single  and  multi-beam  sonar,  mounted  on  a 
vessel.  Relying  on  relatively  high  frequency  sound  to  illuminate  the  bottom,  they  are 
very  accurate  but  time  consuming  owing  to  their  small  swath,  sometimes  taking  many 
days  to  map  a  region  of  interest. 

The  Compact  Hydrographic  Airborne  Rapid  Total  Survey  (CHARTS)  system 
(Figure  1),  an  airplane  with  a  Light  Detection  and  Ranging  (LIDAR)  instrument  mounted 
on  it,  is  also  very  accurate  and  can  perform  large  area  surveys  in  a  much  shorter  time. 
Both  of  these  methods,  however,  require  uncontested  access  to  the  region  of  interest. 
Often  times,  military  operations  are  conducted  in  denied  areas,  where  survey  vessels  and 
aircraft  are  in  danger  of  being  fired  upon.  In  addition,  due  to  the  high  cost  and  demand  of 
these  scarce  resources,  it  is  impractical  to  resurvey  the  same  area  regularly  using  these 
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methods.  This  leads  to  significant  time  intervals  between  surveys  allowing  natural 
nearshore  processes,  such  as  tidal  currents  and  storms,  to  alter  the  bathymetry,  potentially 
rendering  older  surveys  inaccurate. 
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Figure  1.  Compact  Hydrographic  Airborne  Rapid  Total  Survey  (From  Joint 
Airborne  LIDAR  Bathymetry  Technical  Center  of  Expertise,  2011). 
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If  IHO  standards  are  not  required,  other  methods  exist  that  are  better  protected 
from  adversaries  through  greater  stand-off  distances,  are  much  lower  in  cost,  and  have 
frequent  revisit  capability,  especially  those  that  use  satellite-based  sensors.  Nearshore 
bathymetry  can  be  extracted  from  satellite  multispectral  imagery  by  applying  the  linear 
dispersion  theory  of  surface  gravity  waves,  which  is  the  focus  of  this  research.  Previous 
research  into  this  method  at  the  Naval  Postgraduate  School  focused  on  determining  the 
viability  and  accuracy  of  this  method,  but  used  a  manual,  time-intensive  process  (Myrick, 
2011;  McCarthy,  2010).  The  purpose  herein  is  to  take  the  next  step  by  investigating  an 
algorithm  that  automates  this  process.  A  few  of  the  other  methods  that  have  been  or  are 
still  being  explored  are  briefly  discussed  in  Chapter  II. 

B.  SPECIFIC  OBJECTIVES 

The  objective  is  to  investigate  an  automated  method  for  determining  nearshore 
bathymetry  using  remotely  sensed  images  of  the  Waimanalo  Beach  area  of  Hawaii.  A 
desire  to  reduce  the  time  and  effort  required  to  extract  the  bathymetry  infonnation, 
therefore  greatly  increasing  the  efficiency  of  the  process  as  compared  to  the  previous 
work  is  what  provided  motivation  for  this  effort. 

The  WorldView-2  satellite  took  the  multispectral  and  panchromatic  images  in 
rapid  succession  in  March  2011.  The  imagery  was  orthorectified,  converted  to  Universal 
Transverse  Mercator  (UTM)  coordinates,  and  processed  using  a  Wave  Kinematics 
Bathymetry  (WKB)  computer  algorithm  to:  crop  the  images  to  the  desired  analysis  area; 
register  the  images  to  one  another  (~  1  m  accuracy);  apply  filters  and  distinguish  between 
land,  clouds,  and  water  and;  finally,  extract  current  and  depth  fields  by  fitting  a  solution 
to  the  data  (Abileah,  2006).  The  resulting  estimated  depth  fields  were  compared  to 
ground  truth  data  collected  by  the  Joint  Airborne  LIDAR  Bathymetry  Technical  Center  of 
Expertise  (JALBTCX)  during  survey  operations  using  the  CHARTS  system  with  a 
Scanning  Hydrographic  Operational  Airborne  LIDAR  Survey  (SHOALS)-3000 
bathymetric  LIDAR  (figure  1).  This  study  reaffirms  that  applying  the  linear  dispersion 
relation  for  surface  gravity  waves  to  wave  data  collected  from  multispectral  satellite 
imagery  is  a  viable  technique  for  determining  bathymetry  in  denied  or  restricted  areas. 
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II.  BACKGROUND 


A.  INTRODUCTION 

For  many  years,  nearshore  bathymetry  has  been  detennined  without  the  use  of 
direct  measurement,  such  as  a  sounding  line,  even  as  far  back  as  World  War  I  (Myrick, 
2011).  Described  in  more  detail  in  Myrick’s  paper,  those  early  methods  consisted  of  the 
waterline,  transparency,  wave  celerity,  and  wave  period  methods.  The  waterline  method 
relied  on  aerial  photographs  of  the  beach  at  times  of  low  tides,  when  the  bottom  was 
partially  exposed.  The  transparency  method  exploited  the  concept  that  deeper  water 
absorbs  more  light  and,  thus  appears  darker  in  aerial  photographs.  The  wave  period 
method  used  the  fact  that  the  period  is  constant  as  a  wave  propagates  so  a  relationship  can 
be  established  between  the  period,  the  wavelength,  and  the  water  depth  at  multiple 
locations.  These  methods  all  had  limitations  that  made  them  marginally  useful,  but  were 
sometimes  arguably  better  than  having  no  bathymetry  information  at  all  (Myrick,  2011). 

The  wave  celerity  method,  which  invokes  the  linear  dispersion  relation  discussed 
shortly,  was  also  very  limited  before  more  accurate  timekeeping  and  better  resolution 
images  became  available.  Modern  remote  sensing  technology  has  removed  these 
previous  handicaps  and  the  wave  celerity  method  is  now  a  viable  one  for  detennining 
water  depth  in  nearshore  regions  (Myrick,  2011).  The  wave  celerity  method  is  the 
method  used  as  the  foundation  for  this  study. 

B.  THEORY 

Surface  gravity  waves  propagating  in  the  ocean  obey  the  linear  dispersion  relation 
between  wave  celerity,  or  phase  speed,  c;  wave  period  (7);  wavelength  (/.);  and  water 
depth  (d).  The  dispersion  relation  for  surface  gravity  waves  is: 

co1  =  gk  tanh(M) ,  (1) 

where  co  is  the  angular  frequency  (co  =  2k/ I),  g  is  the  acceleration  due  to  gravity,  and  k  is 
the  wavenumber  (k  =  2x/ X)  (Herbers,  2003).  For  large  kd,  which  is  the  case  in  deep 
water,  tanh (kd)  ~  1  and  the  dispersion  relation  reduces  to: 
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CO2  =gk. 


(2) 


In  shallow  water,  where  kd«  1,  tanh (kd)  ~  kd,  so  the  dispersion  relation  reduces  to: 

co-yfgdk .  (3) 

The  corresponding  limits  of  the  phase  speed  c  =  co/k  are: 


CO 

for  deep  water 

(4) 

cs 

II 

toTl 

^4  1 

for  shallow  water. 

(5) 

The  conclusion  from  Equations  (4)  and  (5)  is  that  deep  water  waves  are 
dispersive,  whereas  shallow  water  waves  are  nondispersive.  That  is,  the  phase  speed  of 
deep  water  waves  is  a  function  only  of  co,  whereas  the  phase  speed  of  shallow  water 
waves  is  independent  of  co  and  is  instead,  only  a  function  of  depth  (Herbers,  2003). 
Another  way  to  think  of  this  is,  as  waves  travel  from  deep  to  shallow  water,  their 
frequency  and  period  remain  constant,  forcing  phase  speed  to  decrease  and  wavenumber 
to  increase  with  decreasing  depth.  Since  these  changes  are  proportional  to  each  other, 
this  is  exactly  the  phenomenon  that  is  exploited  to  extract  the  water  depth. 

In  addition  to  waves  slowing  and  getting  shorter  as  they  shoal,  another  physical 
effect  is  an  increase  in  their  amplitude.  This  effect  becomes  important  in  the  surf  zone, 
where  the  waves  break  and  wave  height  increases  wave  speed.  The  nonlinear  processes 
that  occur  in  the  surf  zone  and  add  speed  to  the  waves  are  not  accounted  for  in  the  linear 
dispersion  relation.  This  introduces  error,  usually  in  the  form  of  overestimated  depths 
because  of  the  additional  speed  (Myrick,  2011).  Equation  (1)  is  actually  a  very  good 
approximation  for  depths  greater  than  2  m,  and  is  still  valid  at  1  m  depths  with  moderate 
wave  heights  (Abileah,  2006).  Therefore,  caution  is  necessary  when  using  any  method 
that  relies  on  the  dispersion  relation  in  the  very  shallow  depths  of  the  surf  zone. 

It  can  be  shown  that  the  wave  induced  velocity  components  and  the  water 
pressure  decay  exponentially  with  depth  as  kd  — *  oo  (deep  water).  In  fact,  by  only  half  a 
wavelength  below  the  surface,  these  parameters  are  reduced  to  about  4%  of  their  surface 
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values  (Herbers,  2003).  For  this  reason,  a  good  thumb  rule  is:  regions  where  the  depth  is 
greater  than  half  the  wavelength  are  deep  water  and  the  dispersion  relation  does  not  help 
in  determining  water  depth,  since  phase  speed  does  not  depend  on  depth  in  this  regime. 
For  the  purposes  of  this  study,  regions  where  the  depth  is  less  than  half  the  wavelength 
will  be  considered  “nearshore”  and  is  the  regime  where  the  dispersion  relation  is  useful 
for  this  technique.  It  should  be  noted  that  this  definition  includes  both  the  shallow  water 
regime,  discussed  previously,  and  the  intennediate  water  depth  regime,  where  phase 
speed  depends  on  both  frequency  and  depth. 

C.  SOME  MODERN  METHODS  FOR  DETERMINING  NEARSHORE 

BATHYMETRY 

Several  methods  using  different  techniques  have  been  developed  that  use  modern 
remote  sensing  technology  to  detennine  nearshore  bathymetry.  The  methods  discussed 
are  only  a  few  of  the  many  methods  and  techniques  currently  being  investigated  by  a 
multitude  of  researchers. 

1.  Airborne  Passive  Optical  System  (Using  Linear  Dispersion  Theory) 

The  basis  of  modem  methods  using  the  dispersion  relation  is  to  collect  ocean 
images,  extract  the  space-time  characteristics  of  the  waves  in  the  images,  then  transform 
this  data  into  spectra  that  can  be  used  to  retrieve  depth  and  currents. 

The  first  method  employs  a  passive  optical  system  mounted  on  an  aircraft  (Figure 
2).  This  turret-based  system,  called  the  Airborne  Remote  Optical  Spotlight  System 
(AROSS),  maps  a  time  series  of  images  to  a  common  geodetic  surface  by  carefully 
measuring  the  imaging  geometry  (Dugan  et  al.,  2001a).  It  was  designed  using 
commercial  off-the-shelf  technology  with  the  intention  of  ultimately  mounting  it  on 
unmanned  aerial  vehicles. 

Precisely  registering  the  images  to  each  other  is  critical  to  obtaining  accurate 
results.  Other  important  considerations  are:  adequate  spatial  resolution  to  resolve  the 
smallest  waves  (the  resulting  viewing  geometry  typically  yields  a  sufficient  2  m),  large 
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enough  field  of  view  to  capture  several  of  the  longest  waves  (2  km  x  2  km  is  typically 
used),  and  dwell  times  long  enough  to  observe  several  of  the  longer  wave  periods  (30  s  or 
greater  is  used)  (Dugan  et  al.,  2001a). 


Figure  2.  Close-up  photo  of  camera  and  turret  as  mounted  on  the  nose  cone 
of  the  Pelican  aircraft  (From  Dugan  et  ah,  2001a). 

One  data  set  using  AROSS  was  collected  during  the  Shoaling  Waves  Experiment 
in  the  Outer  Banks  of  North  Carolina  in  1999.  A  three-dimensional  (3-D)  Fourier 
transfonn  was  used  to  turn  a  data  cube  consisting  of  2  minutes’  worth  of  data  from  an 
area  approximately  500  x  500  m  into  a  frequency- wavenumber  power  spectrum.  Figure  3 
shows  a  two-dimensional  (2-D)  slice  through  the  resulting  power  spectrum,  oriented  in 
the  direction  of  the  primary  swell  waves.  The  theoretical  dispersion  surface  that  includes 
no  current  and  is  for  infinite  depth  is  represented  by  the  dashed  curve  at  the  intersection 
of  the  dispersion  surface  with  the  plane  of  the  slice.  The  wave  energy  is  concentrated  in  a 
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narrow  ridge  that  is  slightly  skewed  off  of  the  theoretical  curve.  The  depth  and  current 
are  detennined  by  calculating  a  surface  that  fits  the  measured  wave  spectrum.  The 
resulting  dispersion  surface  is  denoted  by  the  solid  curve  and,  in  this  case,  a  depth  of 
8.2  m  and  current  speed  of  0.14  m/s  is  the  solution  (Dugan,  Piotrowski,  &  Williams, 
2001b). 
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Figure  3.  Frequency- wavenumber  slice,  from  Outer  Banks  study,  through 

the  power  spectrum,  oriented  in  the  direction  of  the  primary  swell  waves. 
The  theoretical  dispersion  surface  that  includes  no  current  and  is  for  infinite 
depth  is  represented  by  the  dashed  curve.  The  depth  and  current  are 
determined  by  calculating  a  dispersion  surface  that  fits  the  measured  wave 
spectrum,  denoted  by  the  solid  curve  (From  Dugan  et  ah,  2001b). 


Another  AROSS  study  was  conducted  in  the  Monterey  Bay  in  1999.  The  same 
size  tiles  were  used,  but  in  this  case,  only  1  minute  of  data  was  Fourier  transformed.  A 
slice  through  the  3-D  spectrum  in  the  direction  of  the  wind  waves  is  shown  in  Figure  4. 
Again,  the  solid  curve  is  the  intersection  of  the  deep  water  dispersion  surface  (i.e.,  no 
current  and  infinite  depth),  but  this  time  the  wave  energy  is  closely  distributed  along  this 
curve.  Since  the  data  closely  match  the  theoretical  curve,  it  can  be  surmised  that  no 
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significant  surface  currents  were  present  during  the  collection.  There  is,  however,  a 
slight  shift  of  the  low  frequency  energy  off  the  curve.  This  is  due  to  the  longer  waves 
(50-100  m  wavelengths)  feeling  the  bottom  because  the  depth  at  the  imaged  location  is 
not  infinite  (i.e.,  is  less  than  half  a  wavelength  deep).  The  actual  depth,  as  detennined 
from  nautical  charts,  is  about  18  m  (Dugan  et  ah,  2001a). 


Wavenumber  (cycles/m) 

Figure  4.  Frequency-wavenumber  slice,  from  Monterey  Bay  study,  through 
the  power  spectrum,  in  the  direction  of  the  wind.  The  solid  curve  is  the 
intersection  of  the  deep  water  dispersion  surface  (i.e.,  no  current  and  infinite 
depth).  The  wave  energy  is  closely  distributed  along  this  curve,  indicating 
that  no  significant  surface  currents  were  present  during  the  collection  (From 

Dugan  et  ah,  2001a). 


One  of  the  conclusions  from  these  and  similar  studies  is  that  a  dwell  (total  image 
time  series  length)  of  1  minute  or  more  is  required  for  high  accuracy.  The  root  mean 
squared  (RMS)  error  approaches  5%  for  2-minute  dwell  times,  but  increases  rapidly  with 
shorter  dwell  times  (Figure  5).  This  graph  is  based  on  using  3-D  Fourier  transforms  in 
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the  inversion  algorithm.  The  algorithm  in  this  study  uses  a  2-D  Fourier  transform  that 
results  in  a  much  flatter  curve  where  errors  do  not  sharply  increase  at  shorter  dwell  times. 
In  the  aforementioned  two  studies,  the  depths  agreed  within  15%  and  the  currents  within 
10%  of  the  in  situ  values  (Dugan,  Piotrowski,  &  Williams,  2002). 
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Figure  5.  RMS  error  as  a  function  of  dwell  (From  Dugan  et  ah,  2002). 

2.  Synthetic  Aperture  Radar  (Using  Wave  Refraction) 

Another  remote  sensing  application  for  determining  nearshore  bathymetry  uses 
Synthetic  Aperture  Radar  (SAR),  which  can  be  done  from  both  airborne  and  space-based 
platforms.  The  basis  for  this  method  is  to  examine  wave  refraction  as  determined  by 
SAR  imagery  using  the  two-scale  Bragg  scattering  model.  The  Bragg  scattering  model 
assumes  that  the  SAR  image  brightness  of  a  patch  of  ocean  is  proportional  to  the 
amplitude  of  the  Bragg  waves.  Bragg  waves  are  small-scale  ocean  surface  waves  that 
have  wavelengths  equal  to  that  of  the  projection  of  the  transmitted  radar  electromagnetic 
wavelength  onto  the  surface,  and  are  propagating  directly  toward  or  away  from  the  sensor 
(Wackerman  &  Clemente-Colon,  2004).  An  expression  for  this  geometry  is: 

Ab=-^—, 

2  sin# 
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where  XB  is  the  Bragg  wavelength,  Xe  is  the  SAR  electromagnetic  wavelength,  and  9  is  the 
local  incident  angle  of  the  ocean  surface.  This  geometry  produces  a  scattering  pattern  of 
the  SAR  radiation  that  leads  to  constructive  interference  and  the  corresponding  brightness 
in  the  image  (Wackerman  &  Clemente-Colon,  2004). 

The  two-scale  aspect  of  the  model  is  illustrated  in  Figure  6,  which  shows  a  slice 
through  a  simplified  ocean  surface  consisting  of  one  large-scale  and  one  small-scale 
wave.  The  Bragg  waves  are  the  smaller  scale  waves  shown  embedded  in  and  riding  on 
the  large-scale  surface  waves.  Flat  plates,  a  few  Bragg  wavelengths  in  size,  are  then  used 
to  model  the  radar’s  interaction  with  the  surface.  The  plates  tilt  and  move  based  on  the 
local  slope  and  motion  of  the  large-scale  wave  surface,  changing  the  radar  image 
signature  as  they  do  (Wackerman  &  Clemente-Colon,  2004). 


Figure  6.  Two-scale  model  illustrating  Bragg  waves  (the  smaller  scale 
waves)  embedded  in  and  riding  on  the  large-scale  surface  waves.  Flat 
plates,  a  few  Bragg  wavelengths  in  size,  model  the  radar’s  interaction  with 
the  surface.  The  plates  tilt  and  move  based  on  the  local  slope  and  motion  of 
the  large-scale  wave  surface  (From  Wackerman  &  Clemente-Colon,  2004). 
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Bathymetry  is  estimated  using  SAR  by  detennining  the  amount  of  wave  refraction 
that  occurs  in  the  images  and  using  propagation  models  that  predict  wave  refraction  as  a 
function  of  depth.  Waves  refract  as  they  move  from  deep  water  toward  the  shore  at  an 
angle  to  the  bathymetric  contours,  which  forces  them  to  turn  and  align  with  the  shore. 
This  turning  is  due  to  a  differential  interaction  with  the  ocean  bottom  (i.e.,  one  end  of  the 
wave  front  reaches  shallower  water  before  the  other  end  and,  therefore,  slows  down 
sooner).  The  faster  the  depth  decreases,  the  more  quickly  the  wave  will  turn  and  become 
parallel  to  the  coast.  Since  refraction  will  not  occur  if  the  waves  start  out  parallel  to  the 
shore  in  deep  water,  this  method  is  ineffective  in  those  instances  (Wackerman  & 
Clemente-Colon,  2004). 


3.  Aerial  Photography  (Using  Radiometric  Techniques) 

A  third  modem  method  relies  on  radiometric  techniques  to  extract  water  depth  in 
shallow  areas  from  multispectral  optical  imagery.  These  techniques  are  based  on 
applying  algorithms  to  interpret  the  radiance  received  by  the  imager.  This  method 
requires  fairly  clear  water  to  work,  so  it  is  not  useful  in  the  many  turbid  nearshore 
locations  around  the  globe. 

A  simple  water  reflectance  model  that  accounts  for  most  of  the  received  signal  is 
represented  by: 


Li  =  LSI  +  krBje 


~Kjfz 


(7) 


where  Lsi  is  the  radiance  over  deep  water;  k,  is  a  constant  that  accounts  for  the  solar 
irradiance,  the  transmittance  of  the  atmosphere  and  the  water  surface,  and  refraction  at 
the  surface;  r«,  is  the  bottom  reflectance;  if  is  the  effective  attenuation  coefficient  for 
water;  /  is  a  geometric  factor  to  account  for  the  path  length  through  the  water;  and  z  is  the 
water  depth  (Lyzenga,  1978).  Using  this  equation  to  solve  for  z  is  the  obvious  way  to 
determine  the  depth. 

A  ratio  algorithm  was  developed  to  help  remove  the  effect  of  changes  in  bottom 
reflectance  on  the  depth  calculation,  so  the  model  can  be  applied  to  scenes  in  which  the 
bottom  composition  changes.  The  algorithm  can  be  represented  by: 
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where  R  is  the  ratio  of  the  bottom-reflected  signals  in  two  bands: 

R  =  (Ll-Lsl)/(L2-Ls2), 


(8) 

(9) 


and  Rb  is  the  ratio  of  the  bottom  reflectances  in  the  same  two  bands.  The  assumption  is 
made  that  the  two  bands  chosen  have  the  same  bottom  reflectance  ratio  for  all  the  bottom 
types  in  the  scene,  which  is  a  decent  assumption  for  multispectral  systems.  This 
algorithm  is  used  with  moderate  success  over  relatively  clear  water  on  data  collected 
from  both  satellites  and  aircraft  (Lyzenga,  1978). 

To  increase  operational  flexibility  and  improve  performance,  a  more  general, 
albeit  more  complex,  algorithm  was  developed  that  modified  the  preceding  model  to 
include  the  effects  of  scattering  in  the  water  and  internal  reflection  at  the  water  surface. 
It  depends  on  proper  classification  of  the  bottom  types  in  the  scene,  which  is  done  prior 
to  the  depth  computation,  and  it  is  at  least  as  good  as  the  ratio  model  at  doing  this  for  the 
case  of  two  bands  (Lyzenga,  1978).  Arguably,  using  more  bands  would  improve  the 
performance  of  the  bottom  classification  portion  of  the  algorithm. 


Figure  7.  Depth  error  comparison  of  the  ratio  and  modified  methods 

(From  Lyzenga,  1978). 
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Once  the  bottom  is  classified,  the  depth,  z,  is  calculated  using: 

Yn  =  Bm  -  Cz ,  (10) 

where  Yn  is  the  set  of  depth  dependent  variables,  Bm  is  a  function  of  the  bottom 
composition  (and  its  value  is  detennined  via  the  aforementioned  classification),  and  C 
characterizes  the  water  attenuation.  This  procedure  introduces  error  both  through  noise 
associated  with  YN  and  if  there  happens  to  be  any  misclassification  of  the  bottom 
(Lyzenga,  1978).  Compared  to  the  ratio  method,  which  is  subject  to  much  higher  noise 
errors,  its  perfonnance  is  much  improved  (Figure  7). 
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III.  EXPERIMENT  DESIGN 


A.  PROBLEM  DEFINITION 

Using  satellite  optical  imagery  to  view  wave  motion  in  order  to  extract  the 
bathymetry  of  nearshore  areas  using  the  linear  dispersion  relation  for  surface  gravity 
waves  has  its  advantages  in  certain  situations.  When  IHO  standards  are  not  required,  it  is 
much  more  cost  effective  and  efficient  to  use  this  method  rather  than  LIDAR-equipped 
aircraft  or  sonar-equipped  vessels.  There  is  no  need  to  dispatch  expensive,  hard-to-get 
assets  to  remote  locations,  since  the  existing  low  earth  orbiting  satellites  will  be  overhead 
most  locations  every  few  days.  It  is  especially  good  for  denied  or  hostile  areas  because 
the  satellites  do  not  attract  attention  and  remain  safely  out  of  the  range  of  most  weapons, 
unlike  aircraft  and  ships. 


Table  1.  Comparison  of  bathymetry  methods  using  remote  sensing 

(From  Abileah,  2006). 


Depth 

accuracy 

Horizontal 

resolution 

Airborne/Satellite 

Turnaround 
(mouths) 1 

Cost  / 
sq  kin1 

Turbid 

waters 

Day/ 

night 

Airborne 

IHO-l 

IHO-l 

Airborne 

6-12 

>$100 

X 

✓ 

Bathymetric  Lidar 

0 

Photobathymetry 

Variable 

10-30  m 

Both 

1-2 

$50 

X 

X 

Currents  interaction 

Variable 

25  nr 

Both 

1 

$10 

✓ 

V 

Wave  refraction 

10-20% 

1  km 

Both 

1 

$10 

✓ 

V 

WKB 

5-10% 

50-200  m 

Both 

1-2 

$50 

✓ 

Notes:  Aough-ordei-of-magntude  assuming  satellites,  except  ABL.  which  is  on  aircraft.  ~Ref:  Argoss  website  (mTO~.argoss.11l  1. 

3WKB  with  the  satellite  imagery  described  in  this  paper  is  daytime  only.  However  nighttime  is  possible  with  long  wave  IR  (Farber. 
1995).  nautical  X-band  radars  (see  text),  and  theoretically  possibly  with  space  SAR. 


Other  satellite  remote  sensing  techniques  offer  these  same  advantages.  However, 
they  each  have  limitations  that  make  each  method  suitable  only  in  certain  situations. 
Table  1  compares  some  of  the  properties  of  each  method.  The  optical  radiance  method 
(Lyzenga,  1978),  referred  to  as  photobathymetry  in  the  table,  is  good  for  relatively  clear 
water  during  cloudless,  daytime  passes  and  has  one  of  the  best  horizontal  resolutions,  but 
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it  is  not  a  good  option  for  turbid  waters,  cloudy  regions,  or  nighttime  passes.  It  also 
requires  in  situ  knowledge  for  calibration  purposes,  negating  its  usefulness  in  inaccessible 
areas  (Abileah,  2006). 

The  wave  refraction  SAR  method  (Wackerman  &  Clemente-Colon,  2004)  is  a 
good  all-weather,  day  or  night  option,  and  it  does  not  rely  on  water  clarity;  but  it  has  the 
least  horizontal  resolution  and  is  not  the  most  accurate  method.  The  WKB  method  in  this 
study,  works  fine  in  turbid  waters,  but  is  not  a  good  option  when  the  scene  is  obscured  by 
clouds  or  for  nighttime  passes  (when  optical  imagery  is  used),  or  when  ocean  waves  are 
not  present  (harbors,  calm  seas,  etc.).  It  has  decent  depth  accuracy  and  horizontal 
resolution,  and  does  not  need  in  situ  data  for  calibrating,  making  it  well  suited  for  denied 
areas.  In  order  to  optimize  the  wave  contrast  in  the  images,  the  camera  should  point  in 
the  direction  of  the  sun  while  avoiding  the  20-30°  sun  glint  cone  (Abileah,  2006). 

This  study  assesses  a  WKB  computer  algorithm  that  automates  the  nearshore 
bathymetry  extraction  on  imagery  of  the  southeastern  portion  of  the  Hawaiian  island  of 
Oahu,  obtained  by  the  WorldView-2  satellite. 

The  S  spectral  bands  of  WorldView-2 

Pan 
Coastal 
Blue 
Green 
Yeltwr 
Red 
Red  Edge 
NIRi 
NIRz 

3  SO  4  SO  550  650  750  £50  S50i  tOSO 

Wave  length  (rim) 

Figure  8.  WorldView-2  sensor  bands  showing  their  relative  positions  and 
overlap  in  the  electromagnetic  spectrum  (From  DigitalGlobe,  2011b). 
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MATERIALS 


1.  WorldView-2  Satellite 

The  WorldView-2  satellite  was  launched  into  a  770  km  high,  sun  synchronous 
orbit  in  October  2009  from  Vandenberg  Air  Force  Base.  The  third  in  DigitalGlobe’s 
commercial  imaging  satellite  constellation,  it  has  a  state-of-the-art  multispectral  optical 
imager  and  provides  approximately  2  m  multispectral  and  0.5  m  panchromatic  resolution 
imagery.  The  panchromatic  and  eight  multispectral  bands  and  their  relative  positions  and 
overlap  in  the  electromagnetic  spectrum  are  illustrated  in  Figure  8.  The  sensor  has  a 
revisit  time  of  1-4  days,  depending  on  desired  viewing  angle  and  resolution.  It  collects 
imagery  in  a  swath  16.4  km  wide  at  nadir  and  can  slew  200  km  in  10  seconds.  Some 
additional  key  specifications  are  listed  in  Table  2  (DigitalGlobe,  2011a). 


Table  2.  Key  WorldView-2  specifications  (After  DigitalGlobe,  2011a) 


Sensor  Bands 

Panchromatic:  450  -  800  nm 

8  Multispectral: 

Coastal:  400  -  450  nm  Red:  630  -690  nm 

Blue:  450  -  510  nm  Red  Edge:  705  -  745  nm 

Green:  510  -  580  nm  Near-IRl:  770  -  895  nm 

Yellow:  585  -  625  nm  Near-IR2  :  860  -  1040  nm 

Sensor  Resolution 

Panchromatic:  0.46  m  GSD  at  nadir,  0.52  m  GSD  at  20°  off-nadir 
Multispectral:  1.85  m  GSD  at  nadir,  2.07  m  GSD  at  20°  off-nadir 

Dynamic  Range 

11-bits  per  pixel 

Swath  Width 

16.4  km  at  nadir 

Attitude  Determination 
and  Control 

3-axis  stabilized 

Actuators:  Control  Moment  Gyros  (CMGs) 

Sensors:  Star  trackers,  solid  state  IRU,  GPS 

Pointing  Accuracy 
and  Knowledge 

Accuracy:  < 500  m  at  image  start  and  stop 

Knowledge:  Supports  geolocation  accuracy  below 

Retargeting  Agility 

Acceleration:  1.43  deg/s/s 

Rate:  3.86  deg/s 

Time  to  Slew  200  km:  10  sec 

Onboard  Storage 

2199  Gb  solid  state  with  EDAC 

Communications 

Image  and  Ancillary  Data:  800  Mbps  X-band 

Housekeeping:  4,  16  or  32  kbps  real-time,  524  kbps  stored,  X-band 
Command:  2  or  64  kbps  S-band 

Max  Viewing  Angle  / 
Accessible  Ground 

Swath 

Nominally  -*-/-45°  off-nadir  =  1651  km  wide  swath 

Higher  angles  selectively  available 

Max  Contiguous 

Area  Collected  in 
a  Single  Pass 
(30®  off-nadir  angle) 

Mono:  138  x  112  km  (8  strips) 

Stereo:  63  x  112  km  (4  pairs) 

Revisit  Frequency 
(at  40*N  Latitude) 

1.1  days  at  1  m  GSD  or  less 

3.7  days  at  20°  off-nadir  or  less  (0.52  m  GSD) 

Geolocation  Accuracy 
(CE.90) 

Specification  of  5.0  m  CE90  at  less  than  30°  off-nadir,  with 
predicted  performance  in  the  range  of  4.6  to  10.7  m  (15 
to  35  ft)  CE90,  excluding  terrain  and  off-nadir  effects 

With  registration  to  GCPs  in  image:  2.0  m  (6.6  ft) 
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2.  SHOALS  Bathymetry 

In  order  to  properly  assess  the  results  obtained  with  the  WKB  algorithm,  a 
comparison  to  ground  truth  depths  was  necessary.  LIDAR  bathymetry,  obtained  for  the 
purpose  of  coral  reef  mapping  by  JALBTCX  using  the  CHARTS  system,  is  provided  as 
ground  truth.  The  survey  was  conducted  in  2000  using  an  Optech,  Inc.,  SHOALS-3000 
LIDAR  instrument  integrated  with  an  Itres  CASI-1500  hyperspectral  imager.  SHOALS 
capabilities  meet  U.S.  Army  Corps  of  Engineers  Hydrographic  Survey  accuracy 
requirements  for  Class  1  surveys  and  the  IHO  nautical  charting  standards  for  Order  1. 
The  positional  accuracy  was  +/-  3  m  in  the  horizontal  and  +/-  15  cm  in  the  vertical  with  a 
resolution  of  0.00000001  degrees  in  latitude  and  longitude  and  0.1  m  in  depth 
(JALBTCX,  2011). 

For  comparison,  it  is  preferable  to  have  recent  data  to  ensure  that  the  bathymetry 
has  not  evolved.  The  SHOALS  data  are  the  most  recent  available  for  the  area  of  interest 
and  was  collected  over  10  years  before  the  satellite  pass  that  produced  the  images 
analyzed  in  this  study.  This  may  be  a  source  of  error,  depending  on  bathymetric 
evolution  during  this  interval.  Another  potential  source  of  comparison  error  is  UTM 
coordinate  in-accuracies,  causing  slight  differences  in  the  pixel  locations  between  the  two 
data  sets.  Another  consideration  is  the  bathymetry  derived  from  the  satellite  data  set 
includes  tide  effects,  whereas  the  SHOALS  data  are  referenced  to  mean  lower  low  water 
(MLLW).  The  satellite  data  must  be  corrected  to  account  for  tide  height.  According  to 
National  Oceanic  and  Atmospheric  Administration  tide  tables,  the  tide  varied  between 
0.06  m  at  low  tide  and  0.50  m  at  high  tide  on  March  31,  2011,  and  was  approximately 
0.37  m  above  MLLW  at  Waimanalo  Beach  at  the  time  of  the  collection. 

3.  WKB  Algorithm 

The  WKB  algorithm  is  a  patent-pending  method  for  generating  maps  of  nearshore 
depth  and  surface  currents  from  a  variety  of  imaging  inputs  from  various  platforms 
(Abileah,  2011).  As  of  December  11,  2011,  it  is  composed  of  almost  100  subroutines 
written  in  MathWorks’  MATLAB  software  and  is  a  work  in  progress  with  several 
updates  planned  over  the  coming  months. 


20 


Figure  9.  STK  snap  shot  of  WorldView-2  collection  pass. 


C.  METHOD 

The  imagery  was  collected  at  about  2200Z  on  March  31,  2011,  by  the 
WorldView2  satellite.  A  Satellite  Tool  Kit  (STK)  snap  shot  shows  the  geometry  part 
way  through  the  collection  (Figure  9).  The  satellite  is  viewing  a  swath  of  the  southeast 
tip  of  Oahu  and  the  adjacent  ocean  from  the  west  looking  toward  the  east  in  the  direction 
of  the  sun,  which  is  the  ideal  geometry  for  high  signal-to-noise  ratio  (SNR)  as  long  as  the 
sun  glint  cone  is  avoided.  The  fact  that  the  waves  are  propagating  east  to  west  toward  the 
sensor  also  improves  SNR.  The  scene  is  mostly  cloud-free  over  the  water  and  has 
adequate  ocean  waves  present,  both  necessary  conditions  for  WKB.  The  Waimanalo 
Beach  area,  covered  by  zone  8  in  the  SFIOALS  bathymetry  (Figure  10),  was  chosen  as 
the  test  area. 

The  imagery  consists  of  two  sets  of  three  pairs  of  panchromatic  and  multispectral 
images  taken  approximately  10  seconds  apart,  and  was  obtained  from  DigitalGlobe  via 
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the  National  Geospatial-Intelligence  Agency.  The  images  were  provided  in  Basic  IB, 
ortho-ready  format  and  geographic  latitude-longitude  coordinates.  An  example  image 
from  the  multispectral  set  is  shown  in  Figure  11.  The  images  were  orthorectified  and 
converted  to  UTM  coordinates  for  the  WKB  algorithm.  A  summary  of  the  multispectral 
images  is  given  in  Table  3.  To  take  advantage  of  their  ability  to  increase  wave  contrast 
by  reducing  subsurface  noise  and  atmospheric  path  radiance,  the  red  and  near-infrared 
bands  (5  through  7)  are  used  for  the  bathymetry  extraction.  The  panchromatic  images, 
which  are  not  listed  in  the  table,  have  the  same  metadata;  however,  the  filenames  contain 
a  “P”  rather  than  an  “M.” 


Figure  10.  SHOALS  2000  survey  of  Oahu  (From  JALBTCX,  2011). 

The  images  are  pre-processed  prior  to  entering  the  WKB  algorithm  by  cropping  to 
a  user-defined  size,  loading  the  cropped  images  into  a  data  cube,  verifying  the  time 
interval  between  successive  images  in  the  set,  registering  the  images  to  each  other, 
conducting  a  deep  water  spectral  analysis,  and  determining  which  pixels  are  water  and 
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which  are  not.  Then  the  program  conducts  a  WKB  extraction  on  the  data  bundle(s)  using 
the  specified  tile  size(s)  (this  determines  the  horizontal  resolution  of  the  resulting  depth 
and  current  fields). 


Table  3.  Summary  of  image  files. 


Image  Name 

First  Line  Time 

Mean  Sat  Az 

Mean  Off-nadir 

View  Angle 

90-mean  Sat  El 

meanSatEl 

Set  1 

11MAR31215059-M1BS-052517305060 

21:50:59 

311.80 

38.5 

44.2 

45.8 

11MAR31215109-M1BS-052517305010 

21:51:09 

306.10 

37.4 

42.8 

47.2 

11MAR31215119-M1BS-0525 17305030 

21:51:19 

300.00 

36.5 

41.8 

48.2 

Set  2 

11MAR31215155-M1BS-052517305050 

21:51:55 

276.30 

35.9 

41.0 

49.0 

11MAR31215204-M1BS-0525 17305040 

21:52:05 

270.10 

36.4 

41.7 

48.3 

11MAR31215214-M1BS-0525 17305020 

21:52:14 

264.20 

37.2 

42.7 

47.3 

Figure  1 1 .  Example  WorldView-2  multispectral  image. 
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1. 


Pre-Processing 


The  pre-processing  ensures  the  images  are  ready  for  the  WKB  algorithm.  The 
current  version  is  described,  but  it  need  not  be  done  exactly  as  it  is  today.  The 
multispectral  and  panchromatic  image  sets  are  processed  separately.  First,  the  image  files 
in  each  set  are  inventoried.  At  least  two  are  required  and  can  be  multispectral  or 
panchromatic  (but  they  have  to  be  the  same  type  with  the  current  version). 

A  graphical  user  interface  (GUI)  allows  several  user  inputs  to  be  entered  during 
the  process.  The  area  for  processing  is  defined  by  cropping  a  Google  Earth  picture 
(Figure  12).  The  cropped  images  are  then  loaded  into  analysis  data  cubes  (one  for  each 
set)  at  a  down-sampled  resolution  that  increases  efficiency  but  still  resolves  the  smallest 
waves  so  as  to  not  lose  any  wave  energy  (3  m  for  the  multispectral  and  1.5  m  for  the 
panchromatic  data).  During  the  loading  process,  the  images  are  stitched  together  to 
create  full  scenes  if  they  were  partitioned  into  smaller  strips  by  the  vendor  (such  was  the 
case  with  this  imagery). 


6.25 
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Figure  12.  Example  cropped  area  (in  green). 
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To  get  the  timing  precisely  correct  for  the  wave  propagation,  the  time  each  scan 
line  was  imaged  is  detennined,  based  on  the  first  line  time  and  the  scan  rate  in  the  image 
metadata.  Without  accurate  timing,  wave  speed  will  not  be  correct  and  the  resulting 
calculations  for  both  current  and  depth  will  be  off. 

Image  registration  is  performed  using  fixed  land  features  to  properly  align  the 
image  frames  with  each  other.  The  area  chosen  should  be  close  to  sea  level  and  have 
buildings  or  other  fixed  structures  in  it  (Figure  13).  Registration  calculates  the  mis¬ 
registration  and  shifts  the  images  to  ensure  the  pixels  at  the  same  geographic  location  on 
each  image  are  overlaid  on  top  of  the  each  other,  to  within  sub-pixel  accuracy.  The 
panchromatic  images  are  used  to  compute  the  registration  shift  because  of  their  higher 
resolution,  but  the  multispectral  could  also  be  used  if  panchromatic  imagery  is  not 
available.  The  multispectral  images  are  shifted  by  the  same  amount  as  the  panchromatic 
ones.  It  should  be  noted  that,  if  the  registration  is  off  slightly,  it  will  simply  lead  to  a 
current  bias  that  does  not  impact  the  bathymetry  detennination. 

Use  mouse  to  outline  registration  chip... 
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Figure  13.  Registration  GUI 
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A  deep  water  box  is  then  analyzed  to  determine  the  time  interval  between  images 
and  to  compare  the  spectra  of  the  multispectral  and  panchromatic  sets  (Figure  14).  The 
interval  is  calculated  using  the  fact  that,  per  the  dispersion  relation  as  applied  to  deep 
water,  wave  speed  is  constant  for  a  given  frequency.  Therefore,  it  is  important  to  pick  a 
spot  where  depth  is  at  least  half  of  a  wavelength  of  the  primary  waves  to  get  an  accurate 
time  interval  (usually  about  50  m).  This  calculation  would  be  unnecessary  if  the  timing 
metadata  were  more  precise.  It  is  important  to  determine  the  timing  very  precisely  for  the 
wave  speed  calculation.  The  interval  timing  is  combined  with  the  time  offset  determined 
earlier  to  obtain  very  accurate  time  spacing  for  pixels  in  different  images.  Other  analyses 
are  conducted  with  regard  to  radiance  that  are  not  important  for  this  work,  but  may  be  of 
interest  to  researchers  using  radiometric  techniques. 
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Figure  14.  Deep  water  analysis 


The  final  pre-processing  function  creates  a  mask  that  prevents  pixels  that  are  not 
water  from  being  included  in  the  spectral  analysis  (Figure  15).  The  various  bands  of  the 
multispectral  images  are  compared  to  determine  if  a  pixel  most  likely  contains  land, 
water,  clouds,  a  boat  or  some  other  interfering  phenomenon.  With  the  “bad”  pixels 
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removed,  and  all  the  other  adjustments  made  during  pre-processing,  the  images  are  ready 
to  go  into  the  WKB  algorithm  for  extraction  of  the  bathymetry. 

2.  WKB  Extraction 

WKB  is  rooted  in  the  fact  that  surface  gravity  waves  decrease  in  speed  as  they 
move  to  shallower  depths,  according  to  the  linear  dispersion  relation  given  by  Equation  1 . 
By  capturing  images  of  ocean  waves  at  short,  precisely  known  time  intervals,  their 
horizontal  displacement  can  be  observed  and  celerity  calculated.  Depth  is  inferred  from 
the  dispersion  relation.  WKB  becomes  less  useful  in  the  cases  where  depths  approach  the 
height  of  the  waves,  such  as  in  the  surf  zone  (due  to  non-linear  effects  that  are  not 
incorporated  in  the  linear  dispersion  relation);  and  where  depths  are  greater  than  half  the 
wavelength  of  the  longest  detectable  waves  (the  region  where  wave  speed  is  independent 
of  depth). 


Figure  15.  True  color  image  (left)  with  corresponding  land  mask  (right) 

The  key  to  extracting  bathymetry  from  the  information  contained  in  the  images  is 
the  Fourier  transform.  Most  commonly,  3-D  Fourier  analysis  is  applied  to  transform 
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image  data  from  two-dimensional  space-time  into  wavenumber-frequency  spectra.  This 
technique  requires  long  time  series  lengths,  or  dwell  times  -100  seconds  (Dugan  et  ah, 
2001a). 

Two-dimensional  Fourier  analysis  is  more  appropriate  for  satellite  imagery, 
which,  due  to  the  rapidly  changing  geometry,  has  effective  dwell  times  equal  to  the  time 
between  images.  The  2-D  analysis  allows  for  these  much  shorter  ~10  second  dwell  times 
without  adversely  increasing  the  error.  This  is  because  the  depth  accuracy  does  not 
depend  on  dwell  time  so  the  frequency  resolution  is  not  relevant  (Abileah  &  Trizna, 

2010) .  In  this  case,  the  2-D  data  are  transformed  into  2-D  wavenumber  spectra.  In  both 
methods,  the  image  is  broken  into  2-D  tiles,  the  dimensions  of  which  are  driven  by  a 
compromise  between  maximizing  spatial  resolution  and  minimizing  depth  error  (Abileah, 

2011) .  Tiles  that  are  too  small  provide  better  resolution  at  the  expense  of  depth  accuracy. 
Large  tiles  produce  better  accuracy,  but  resolution  suffers  since  one  depth  is  calculated 
per  tile. 

The  2-D  algorithm,  a  modification  of  which  is  applied  in  this  study,  applies  a 
propagation  kernel  to  the  2-D  Fourier  transforms  of  the  N  images  {^(x^S^x),...,  SN  (x)} 

,Fn(k)  -  32[,SJ  .  The  propagation  kernel  moves  the  waves  forward  or  backward  in  time 
and  is  defined  as, 


0±(kx,ky\d,u)  =  e 


±ico0(kx,ky\d,u)T 


(11) 


where  d  is  depth,  u  is  the  surface  current,  r  is  the  time  interval  between  compared 
images,  and  kx  and  ky  are  the  x  and  y  components  of  the  wavenumber.  The  sign  in  the 
exponent  detennines  whether  the  waves  propagate  forward  or  backward.  The  term 
co0(kx,kv  |  d,u)  in  the  exponent  comes  from  the  linear  gravity  wave  dispersion  relation  in 
this  form: 


co0(k  \d,u)  =  |  k  |  tanh(|  k\d)-k  ■ u 


(12) 
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which  also  accounts  for  modification  of  the  waves  by  ocean  currents.  Using  the 
propagation  kernel,  the  Fourier  transfonn  at  time  n  is  related  to  the  next  at  time  n+m  by 

Fn+m  =  0"'F  (Abileah  &  Trizna,  2010).  (13) 

Wave  energy  in  the  transform  space  tends  to  be  distributed  along  a  dispersion 
surface  and  the  depth  and  current  can  be  found  by  finding  the  best  fit  of  the  dispersion 
formula  to  this  measured  power  spectrum.  One  way  to  find  the  best  fit  is  by  minimizing 
the  difference  between  successive  Fourier  transfonns  by  tuning  the  propagation  kernel  to 
the  correct  depth  and  current.  Mathematically,  this  is  done  by  minimizing  the  objective 
function 

N—m 

argmin,  ,  JA  =  £  W(k)^  I  F„„  |2 

,  (14) 

where  W ( k )  is  a  weighting  filter  that  increases  the  accuracy  of  the  solution  by  increasing 
the  SNR.  It  does  this  by  masking  much  of  the  background  noise  by  keeping  those 
wavenumbers  where  the  wave  energy  is  concentrated  and  eliminating  those  where  it  is 
not.  Since  waves  are  located  where  the  energy  adds  coherently  in  wavenumber  space 
when  the  images  are  propagated  to  the  same  time,  this  determination  is  done  by  seeing 
for  which  wavenumbers  this  is  true.  The  objective  function  also  includes  a  provision  for 
more  than  two  image  pairs  by  including  a  summation  over  frame  intervals,  in,  from  n  =  1 
to  N-m  (Abileah,  2011). 

A  patent-pending  technique,  which  is  a  modification  of  the  2-D  algorithm, 
overcomes  the  accuracy-resolution  compromise  discussed  above  by  adding  a  2-D  inverse 
Fourier  transform  operation,  3  j1 ,  to  the  objective  function  (Abileah,  2011).  By 

transfonning  the  data  from  wavenumber  space  back  into  spatial  coordinates,  a  depth  is 
obtained  for  every  image  pixel  rather  than  for  every  tile.  In  fact,  it  eliminates  the  need  to 
tile  the  images,  although  this  has  not  yet  been  incorporated  into  the  code.  To  expedite 
implementing  the  new  method,  the  existing  code  was  modified  rather  than  generating 
new  code,  which  is  computationally  inefficient  but  yields  the  same  results.  A  more 
efficient  revision  is  forthcoming. 


29 


Besides  better  spatial  resolution  and  less  computation  time,  another  benefit  of  the 
new  method  is  it  is  less  susceptible  to  interference  from  non-water  pixels  because,  unlike 
in  wavenumber  space,  in  spatial  coordinates  the  waves  can  be  separated  spatially  from 
the  interfering  phenomena. 

As  before,  the  2-D  Fourier  transfonn  is  applied  to  the  (in  this  case,  tiled)  images, 
converting  them  into  2-D  wavenumber  spectra.  (New  implementations  can  forgo  the 
tiling  step;  however,  tiling  may  still  be  desired  to  take  advantage  of  the  multiple 
processors  in  parallel  computing  settings.)  Minimizing  the  revised  objective  function: 


ar§  Ja  0)  =  32!  \_w(k )(Fi  (k )  -  °  'F2 (k )) 


(15) 


again  by  tuning  ®,  produces  the  best  depth  and  current  solution  for  each  x  (Abileah, 

2011). 

Although  W  eliminates  much  of  the  noise,  some  is  still  present  after  the  filtering. 
To  further  improve  the  results,  a  sliding  window  A(x )  is  applied.  The  window  size  can 
be  adjusted  as  necessary  to  improve  SNR;  however,  a  larger  window  will  reduce 
resolution,  but  not  to  the  degree  the  tiles  did  in  the  old  method.  With  the  appropriate 
A(x) ,  spatial  resolution  with  the  new  method  can  be  10  m  (Abileah,  2011).  The  code 
implemented  herein  used  250  m  tiles  with  50%  overlap,  creating  an  effective  spatial 
resolution  of  about  125  m. 

For  instances  where  more  than  two  images  are  available,  SNR  can  be  increased 
further  by  averaging  across  the  images,  n  =  1  to  N-m,  not  unlike  in  the  old  method.  For 
data  sets  with  long  interval  times,  and  as  was  done  for  this  study  with  10  seconds 
between  images,  consecutive  images  are  used,  meaning  m=  1 .  Still  another  improvement 
through  averaging  is  to  sum  over  a  set  of  interval  spacing  values,  M,  so  multiple  image 
separation  options  are  available.  The  final,  general  form  for  the  objective  function  that 
includes  all  these  possibilities  is  (Abileah,  2011): 

argmin dJlJA(x)=  £  £  £  [W{k){F, (k)-®-\F2(k))]|  .  (16) 

me  A/  n= 1  A(x) 
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Choosing  an  m  other  than  one  creates  comparisons  of  images  that  are  not 
consecutive.  Since  this  increases  the  time  interval  between  compared  images,  effectively 
decreasing  the  sampling  rate,  it  should  be  done  with  the  following  consideration  in  mind. 
In  wavenumber  space,  a  phase  shift  is  equivalent  to  wave  displacement  in  spatial 
coordinates.  Using  Equation  11,  the  measured  phase  shift  leads  to  the  depth  result 
obtained  with  Equation  16.  In  theory,  a  very  small  time  interval  is  all  that  is  necessary  to 
determine  the  phase  shift.  In  the  real  ocean,  where  noise  complicates  the  wave  signal, 
however,  a  longer  time  is  required  to  accurately  measure  the  phase  shift  amidst  the 
background  noise.  On  the  other  hand,  too  long  of  an  interval  allows  the  waves  to  become 
less  coherent  as  they  become  modified  by  wind,  refraction,  bottom  friction,  etc. 
Therefore,  a  time  interval  of  at  least  one  second,  but  not  more  than  several  seconds,  is  a 
good  compromise.  Ideally,  the  10-second  interval  time  between  the  images  in  this  study 
should  be  shorter  to  ensure  less  wave  modification.  With  this  in  mind,  an  m  larger  than 
one  would  only  be  appropriate  for  data  sets  having  image  interval  times  (or  sampling 
times)  of  one  second  or  less  between  images.  This  can  also  be  extended  to  using  sets  of 
different  values  of  m  (Abileah,  2011). 

There  are  a  couple  of  ways  to  implement  the  objective  function  tuning,  which 
involves  searching  for  the  depths  and  current  speeds  that  minimize  Equation  16.  One 
method  attempts  this  by  trying  every  possible  combination  of  depth  and  current  speed. 
This  option  is  computationally  intensive  since  it  has  to  loop  through  three  variables, 
depth  and  both  the  x  and  y  components  of  the  current  velocity,  in  computing  the 
propagation  kernels  for  use  by  Equation  16  (Abileah,  2011). 

Another  way  to  do  it  that  reduces  the  computation  time  significantly  (and  the  way 
the  current  revision  accomplishes  it)  is  to  separate  the  current  search  from  the  depth 
search.  To  do  this,  a  low  wavenumber  filter  is  applied  to  isolate  the  higher  wavenumbers 
(short  wind  waves)  because  they  are  effectively  deep  water  waves  and  unaffected  by 
depth.  This  removes  the  need  to  include  depth  in  the  first  search  and  allows  any 
deviation  from  the  dispersion  relation  predicted  wave  propagation  to  be  attributed  to  the 
current.  Once  the  current  is  known,  the  images  are  shifted,  similar  to  the  registration 
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process,  to  remove  any  wave  displacement  caused  by  the  current.  Then  the  search  for 
depth  is  done  to  complete  the  solution  (Abileah,  2011). 

Of  note,  in  execution  of  Equation  13,  the  current  code  does  not  propagate  the 
waves  in  one  image  to  the  time  of  the  other  image  being  compared.  Rather,  it  propagates 
the  waves  in  both  images,  one  forward  and  one  backward,  to  the  time  half-way  between 
the  two  (Abileah,  2011).  This  was  done  purely  for  cosmetic  reasons. 

The  output  of  the  WKB  program  includes  several  figures,  some  of  which  are 
shown  above,  a  diary  that  contains  a  chronological  account  of  the  program’s  operations, 
image  registration  movies  that  show  the  shift  being  applied  to  the  images,  and  image  data 
cubes  that  contain  the  image  data  at  various  stages  of  the  process.  The  main  output  figure 
is  a  quad  chart  that  displays  wave  direction  on  the  top  left,  a  true  color  multispectral 
image  of  the  area  processed  on  the  top  right,  the  extracted  bathymetry  on  the  bottom  left, 
and  the  extracted  ocean  currents  on  the  bottom  right  (Figure  16).  Another  key  figure, 
which  will  be  discussed  in  Chapter  IV,  shows  how  the  WKB  derived  bathymetry 
compares  to  the  SHOALS  bathymetry. 


32 


UTM  Northing  (km)  UTM  Northing  (km) 


Peak  wave  direction  &  log  wavelength 


Image 


2364 


2363 


2362 


2361 


2360 


2359 


2358 

634  635  636  637  638  639  640 

UTM  Easting  (km) 


2364 


■g  2363 
<?2362 
I  2361 

Z 

f!  2360 
2359 


2358 


634  635  636  637  638  639  640 

UTM  Easting  (km) 


Depth 


2364 


2363 


2362 


2361 


2360 


2359 


2358 

634  635 


636  637  638  639  640 

UTM  Easting  (km) 


Currents 


2364 

2363 

2362 

2361 

2360 

2359 

2358 


634  635  636  637  638  639  640 


UTM  Easting  (km) 


Figure  16.  Image  Set  1  case  3MSI  WKB  output  showing  wave  direction  (top 
left),  a  true  color  image  (top  right),  extracted  bathymetry  (bottom  left),  and 
extracted  ocean  currents  (bottom  right) 
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IV.  ANALYSIS  RESULTS 


A.  GENERAL  APPROACH 

For  the  comparison  analysis  of  WKB  extracted  bathymetry  to  SHOALS 
bathymetry,  two  data  runs  were  conducted  by  ingesting  two  different  image  sets  into  the 
WKB  algorithm.  The  first  set  consisted  of  the  first  three  multispectral  images  in  Table  3 
and  their  corresponding  panchromatic  images,  while  the  second  set  consisted  of  the  last 
three  multispectral  images  and  their  corresponding  panchromatic  images.  While  it  was 
not  practical  to  make  the  user-selectable  items  identical  between  runs,  the  cropped  area  of 
interest,  registration  box,  and  deep  water  box  were  carefully  chosen  to  minimize 
differences  between  the  two  sets.  This  way,  any  differences  in  the  results  would  be 
mostly  attributable  to  differences  in  the  images. 

Four  cases  were  produced  in  each  of  the  two  data  runs,  two  using  multispectral 
images  and  two  using  panchromatic  images.  The  four  cases,  named  for  how  many  of 
each  type  of  image  were  ingested,  are:  two  multispectral  images  (2MS1),  three 
multispectral  images  (3MS1),  two  panchromatic  images  (2Pan),  and  three  panchromatic 
images  (3Pan).  For  the  two-image  cases,  the  first  two  images  in  each  set  were  the  ones 
ingested. 

The  WKB  extracted  depth  fields  were  compared  to  the  SHOALS  bathymetry 
using  side-by-side  area  maps.  These  figures  show  the  SHOALS  map  on  the  left  and  the 
WKB  map  on  the  right  and  allow  a  quick  visual  comparison  of  the  WKB  field  with 
ground  truth  field  (Figure  17). 

East-west  transects  were  chosen  to  investigate  a  few  sample  cross-sections  for 
each  case.  For  each  transect,  the  WKB  depth  is  plotted  with  the  SHOALS  depth  as  a 
measure  of  the  similarity  between  them  along  the  east- west  slice.  The  depth  error  is 
calculated  for  each  transect,  as  well.  This  graph  aids  in  identifying  where  along  the 
transect  the  greatest  deltas  are  located.  A  third  graph  for  each  transect  is  the  relative 
depth  error,  which  shows  depth  error  as  a  percentage  of  depth.  This  allows  visualization 
of  where  the  error  has  more  or  less  significance  (Figure  18). 
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Figure  17.  Comparison  maps  of  SHOALS  (left)  and  3MSI  WKB  for  Image 

Set  1  (right) 


In  addition  to  transects,  each  case  was  analyzed  using  all  the  valid  depths  in  the 
entire  subject  area.  Scatter  plots  show  the  distribution  of  the  data  and  their  correlation. 
Since  the  data  field  consists  of  millions  of  points,  resulting  in  saturated  plots,  the  data 
were  thinned  out  by  picking  out  every  thousandth  point  to  plot.  The  plots  specify  the 
slope  of  the  linear  regression  line,  its  y-intercept,  and  R  ,  the  square  of  the  correlation 
coefficient  for  the  SHOALS-observed  and  WKB-predicted  values. 

The  regression  line  slope  and  R2  are  indicators  of  WKB  performance.  If  WKB 
were  perfect,  its  depths  would  be  exactly  correlated  to  the  SHOALS  depths,  which  would 
yield  slope  and  R  values  of  1.0.  The  less  perfectly  WKB  performs,  the  less  correlated 
the  WKB  and  SHOALS  depths  will  be,  and  the  further  these  values  decrease  from  1.0. 
The  scatter  plots  have  a  second  panel  that  filters  out  all  SHOALS  depths  greater  than  15 
m,  for  reasons  discussed  later  (Figure  19). 
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Figure  18.  Transect  from  Image  Set  1  case  3MSI  showing  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


Figure  19.  Scatter  plot  from  Image  Set  1  case  3MSI  showing  thinned  data  for 
all  depths  (top)  and  for  just  depths  less  than  15  m  (bottom) 
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In  addition,  a  bar  graph  was  produced  that  shows,  for  each  case,  the  mean  error 
for  several  depth  partitions,  or  bins.  The  bins  are  one  to  2  m,  followed  by  every  2  m  up 
to  20  m,  with  the  last  bin  including  all  depths  greater  than  20  m.  This  aids  in  determining 
at  which  depths  the  algorithm  worked  better  (Figure  20). 


WKB  Depth  Error  vs  Depth 


Depth  (m) 

Figure  20.  Bar  graph  from  Image  Set  1  case  3MSI  showing  the  mean  depth 

error  for  several  depth  bins 


B.  IMAGE  SET  1 

A  quick  glance  at  the  comparison  map  for  3  MSI  (Figure  17)  shows  that,  overall, 
WKB  did  a  good  job  of  depicting  the  main  bathymetric  features,  such  as  the  elongated 
trough  (yellow  feature)  near  the  center  of  the  map  and  the  shallower  spots  to  the  south  of 
it  (dark  red).  The  abrupt,  steeper  slope  beyond  the  more  gently  sloping  beach  is  also  well 
depicted  where  the  red  quickly  turns  yellow  a  couple  of  kilometers  from  the  beach. 
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Image  Set  1  results  are  summarized  in  a  line  graph  that  consolidates  the  mean 
depth  errors  in  the  bar  graphs  of  all  four  cases  on  one  chart  (Figure  21),  in  a  table  listing 
the  mean  depth  error  and  relative  error  ranges  for  each  depth  bin  (Table  4),  and  in  a  table 
listing  the  R2  values  for  each  case  (Table  5).  The  results  highlight  some  of  the  strengths 
and  weaknesses  of  the  method.  The  complete  series  of  graphs  for  Image  Set  1  is  found  in 
Appendix  A. 


WKB  Depth  Error  vs  Depth 


Figure  2 1 .  Mean  depth  errors  for  all  four  Image  Set  1  cases 

Understandably,  the  1-2  m  bin  for  all  cases  is  dominated  by  extremely  large  mean 
depth  errors.  The  nonlinear  processes  in  the  surf  zone  render  the  linear  dispersion 
relation  inaccurate,  and  therefore,  WKB  less  effective  at  predicting  depths  below  about 
2  m.  To  improve  WKB  in  this  regime  would  require  it  to  account  for  nonlinearities  as  the 
depth  approaches  the  height  of  the  waves.  This  could  be  done  by  applying  a  correction 
within  the  linear  dispersion  relation  (Fledges,  1976),  but  would  still  be  an  approximation. 
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Other  factors  contributed  to  error  at  these  shallow  depths,  including  an  ambiguity  caused 
by  a  sampling  time  that  is  equivalent  to  the  period  of  some  of  the  waves,  and  insufficient 
short  waves  for  current  estimation. 


Table  4.  Mean  depth  error  and  relative  error  ranges  for  each  depth  bin  in  Image  Set  1 


Case 

<2m 

2-4m 

4- 6m 

6-8m 

2MSI 

6.8 

341% 

3.8 

94% 

188% 

0.9 

14% 

21% 

1.0 

13%  17% 

3MSI 

8.4 

420% 

1.6 

40% 

80% 

0.6 

11% 

16% 

0.8 

10%  14% 

2Pan 

3.9 

196% 

1.7 

42% 

84% 

1.3 

22% 

33% 

1.9 

24%  32% 

3Pan 

5.8 

288% 

1.1 

28% 

55% 

1.2 

20% 

29% 

1.3 

16%  22% 

Case 

8- 10m 

10-12m 

12- 14m 

2MSI 

1.4 

14% 

17% 

2.3 

19% 

23% 

3.6 

26% 

30% 

3MSI 

1.5 

15% 

18% 

2.5 

21% 

25% 

2.1 

15% 

18% 

2Pan 

3.0 

30% 

37% 

2.4 

20% 

24% 

3.4 

24% 

29% 

3Pan 

2.9 

29% 

37% 

3.1 

26% 

31% 

1.9 

14% 

16% 

Case 

14- 16m 

16- 18m 

18-20m 

>20m 

2MSI 

4.7 

29% 

33% 

4.5 

25% 

28% 

5.6 

28% 

31% 

4.7 

23% 

3MSI 

3.8 

24% 

27% 

4.3 

24% 

27% 

5.6 

28% 

31% 

4.8 

24% 

2Pan 

4.6 

29% 

33% 

6.6 

37% 

41% 

6.6 

33% 

37% 

5.2 

26% 

3Pan 

3.5 

22% 

25% 

4.8 

27% 

30% 

6.1 

31% 

34% 

4.8 

24% 

2 

Table  5.  R~  values  for  Image  Set  1 


Case 

R2 

R2  <  15  m 

2MSI 

0.8681 

0.8909 

3MSI 

0.9005 

0.9078 

2Pan 

0.8603 

0.8453 

3Pan 

0.905 

0.911 
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Another  limitation  occurs  at  the  other  end  of  the  depth  field,  where  mean  depth 
errors  were  relatively  high.  At  depths  greater  than  half  the  wavelength,  waves  are  not 
affected  significantly  by  the  bottom,  so  they  are  dispersive  and  the  dispersion  relation  is 
not  helpful  for  determining  depth.  Therefore,  WKB  can  only  determine  depth  up  to  the 
point  where  the  water  depth  allows  gravity  waves  to  become  nondispersive.  If  fewer 
longer  waves  are  present  in  the  wave  field,  this  will  occur  at  shallower  depths.  If  longer 
waves  are  abundant,  then  WKB  will  be  effective  at  deeper  depths. 

Therefore,  WKB  is  limited  on  its  ability  to  extract  depth  as  determined  by  the 
wave  spectrum  detected  by  the  sensor.  In  this  instance,  it  appears  that  the  spectrum  lacks 
sufficient  energy  in  waves  longer  than  about  40  m,  so  that  by  a  depth  of  about  16  m,  the 
mean  depth  errors  are  consistently  above  4  m  with  relative  errors  of  25%  and  greater. 
This  assumption  is  supported  by  The  Coastal  Data  fnfonnation  Program  (CD1P)  of 
Scripps  Institution  of  Oceanography,  which  indicates  the  sea  on  3 1  March  off  of  nearby 
Mokapu  Point  was  relatively  calm  with  very  little  long  wave  swells  (Figures  22  and  23). 
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Figure  22.  Energy  spectrum  from  Mokapu  Point  buoy  for  March  2011, 
showing  larger  swell  events  earlier  in  the  month,  but  relatively  little  swell  on 
the  date  of  collection  (31  March)  (From  CDIP,  2012) 


42 


Figure  23.  Energy  spectra  from  Mokapu  Point  buoy  for  the  time  of  collection 
(left)  and  a  swell  event  a  few  days  prior  (right),  showing  an  order  of 
magnitude  difference  in  energy  at  the  low  frequencies  (From  CDIP,  2012) 

In  between  these  limits,  perfonnance  was  much  better.  3MSI  had  the  best  overall 
performance,  with  a  mean  depth  error  of  1.8  m  from  2-16  m,  a  regression  line  slope  of 
1.01  and  R2  of  91%.  The  next  best  was  3Pan  with  2.1  m  mean  depth  error,  slope  of  1.06 
and  R "  of  91%.  The  two-image  cases  did  the  worst  with  mean  depth  errors  of  2.5  m  and 
2.6  m  for  2MSI  and  2Pan,  respectively.  It  was  expected  that  the  three-image  cases  would 
perform  better  than  the  two-image  ones  because  of  the  averaging  technique  that  increases 
SNR  with  multiple  image  pairs.  However,  the  multispectral  outperfonning  the 
panchromatic  was  unexpected  and  is  addressed  shortly. 

In  the  shallower  depths  outside  the  surf  zone,  the  multispectral  cases  performed 
quite  well.  From  about  four  to  10  m,  mean  depth  errors  were  1.5  m,  or  less,  and  relative 
errors  were  in  the  10-20%  range.  By  10  m,  the  best  panchromatic  case  (2Pan  in  this 
depth  bin)  did  as  well  as  the  multispectral  cases,  all  having  mean  depth  errors  of  about 
2.4  m  and  relative  errors  in  the  20-25%  range.  A  bit  of  an  anomaly,  this  was  the  only 
instance  where  2Pan  was  better  than  3Pan.  Looking  at  the  comparison  maps  in  Appendix 
A,  it  may  be  due  to  3Pan  doing  a  slightly  worse  job  in  the  area  west  of  Manana  Island  in 
the  southeast  corner,  an  area  which  is  primarily  10-12  m  in  depth. 
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The  expectation  was  that  the  panchromatic  cases  would  outperfonn  the 
multispectral  cases  because  of  their  higher  resolution  images  being  able  to  better  resolve 
the  shorter  waves  leading  to  a  more  accurate  current  determination.  This,  in  turn,  would 
produce  more  accurate  depths.  A  possible  reason  for  the  multispectral  cases’  better 
performance  has  to  do  with  water  clarity.  The  ocean  bottom  can  be  seen  quite  clearly  in 
the  shallower  regions  in  Figure  11.  This  certainly  added  to  the  non- gravity  wave  noise  in 
the  panchromatic  cases  with  a  corresponding  reduction  in  SNR.  The  multispectral  cases 
were  partially  immune  to  this  problem  by  screening  out  most  of  the  subsurface  noise  by 
using  only  bands  5  through  7.  What  the  multispectral  cases  lacked  in  resolution,  they 
made  up  for  in  band  selectivity  and  vice  versa  for  the  panchromatic  cases.  In  the  deeper 
water,  where  bottom  radiance  was  less  significant,  the  panchromatic  cases  did  better. 

In  line  with  this  reasoning,  3Pan  perfonned  the  best  in  the  12-16  m  region,  edging 
out  3MSI  by  a  couple  tenths  of  a  meter  in  depth  error.  Mean  depth  errors  were  about  2  m 
in  the  12-14  m  bin  and  3.5  m  in  the  14-16  m  bin  for  these  two  cases.  The  worse¬ 
performing  two-image  cases  had  well  over  3  m  errors  in  the  12-14  m  bin  and  well  over  4 
m  in  the  14-16  m  bin. 

C.  IMAGE  SET  2 

Image  Set  2  yielded  similar  results  to  set  one,  with  some  notable  differences. 
Again,  WKB  did  a  good  job  of  depicting  the  main  bathymetric  features,  as  seen  in  the 
comparison  maps  in  Appendix  B,  where  a  complete  series  of  Image  Set  2  graphs  are 
found.  As  was  done  for  Image  Set  1,  Image  Set  2  results  are  summarized  in  a  line  graph 
of  all  four  cases  that  consolidates  the  mean  depth  errors  on  one  chart  (Figure  24),  in  a 
table  listing  the  mean  depth  error  and  relative  error  ranges  for  each  depth  bin  (Table  6), 
and  in  a  table  listing  the  R~  values  for  each  case  (Table  7). 

The  same  limitation  in  the  surf  zone  is  apparent  in  Image  Set  2’s  results,  with  the 
1-2  m  bin  for  all  cases  dominated  by  extremely  large  mean  depth  errors.  In  the  deeper 
region,  although  the  same  limitations  applied  to  this  set,  better  results  were  achieved  with 
set  two.  For  instance,  a  respectable  mean  depth  error  of  3.0  m  and  corresponding  relative 
error  of  17-19%  was  obtained  by  3MSI  in  the  16-18  m  bin.  In  comparison,  all  cases  in 
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set  one  had  errors  well  above  4.0  m  and  greater  than  23%  by  16  m  depth.  In  fact,  3MSI 
performed  better  in  the  18-20  m  bin  in  set  two  than  it  did  in  the  16-18  m  bin  in  set  one. 
It  was  not  until  the  final  greater-than-20  m  bin  that  3MSI’s  performance  degraded  above 
the  poorer  set  one  values.  The  performance  may  have  been  aided  by  geometry  more 
favorable  for  enhancing  long  wave  contrast.  If  this  was  the  case,  the  wave  spectrum  has  a 
higher  SNR  in  the  long  wave  portion,  allowing  for  deeper  bathymetry  extraction. 
However,  this  is  purely  speculative. 


WKB  Depth  Error  vs  Depth 


Figure  24.  Mean  depth  errors  for  all  four  Image  Set  2  cases 

As  in  set  one,  the  three-image  cases  outperformed  their  two-image  counterparts 

almost  exclusively.  Also  like  set  one,  3MSI  had  the  best  overall  performance  in  between 

the  surf  zone  and  the  deeper  limit,  with  a  mean  depth  error  of  1.7  m  from  2-16  m,  and  an 

only  slightly  larger  1.9  m  when  expanded  to  include  16-18  m.  Regression  line  slope  and 

R~  for  3MSI  were  1.13  and  88%,  respectively.  The  next  best  for  set  two  was  2MSI, 

rather  than  3Pan  as  was  the  case  with  Image  Set  1 .  It  had  a  mean  depth  error  of  2.0  m  for 
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the  range  of  2-16  m  and  similar  regression  line  slope  and  R  .  The  panchromatic  cases  did 
the  worst  this  time  with  2.6  m  and  3.4  m  for  3Pan  and  2Pan,  respectively.  These  figures 
are  quite  a  bit  worse  than  the  2.1m  and  2.6  m  achieved  with  set  one. 


Table  6.  Mean  depth  error  and  relative  error  ranges  for  each  depth  bin  in  Image  Set  2 


Case 

<2m 

2-4m 

4- 6m 

6-8m 

2MSI 

7.4 

369% 

1.6 

40% 

80% 

0.8 

13% 

20% 

1.0 

13% 

17% 

3MSI 

5.1 

256% 

1.6 

41% 

82% 

0.7 

11% 

16% 

1.4 

18% 

23% 

2Pan 

18.1 

903% 

4.4 

111% 

222% 

2.8 

47% 

71% 

3.4 

43% 

57% 

3Pan 

14.5 

726% 

3.1 

77% 

155% 

3.2 

54% 

80% 

2.7 

34% 

45% 

Case 

8- 10m 

10- 12m 

12- 14m 

2MSI 

2.6 

26% 

33% 

2.3 

19% 

23% 

2.0 

14% 

17% 

3MSI 

2.2 

22% 

27% 

1.8 

15% 

18% 

1.8 

13% 

15% 

2Pan 

3.8 

38% 

48% 

2.9 

24% 

29% 

3.5 

25% 

29% 

3Pan 

3.6 

36% 

45% 

1.7 

14% 

17% 

1.5 

11% 

13% 

Case 

14- 16m 

16- 18m 

18-20m 

>20m 

2MSI 

3.4 

21% 

24% 

3.8 

21% 

24% 

4.6 

23% 

25% 

5.4 

27% 

3MSI 

2.4 

15% 

17% 

3.0 

17% 

19% 

3.9 

20% 

22% 

5.0 

25% 

2Pan 

3.2 

20% 

23% 

5.3 

29% 

33% 

5.4 

27% 

30% 

5.0 

25% 

3Pan 

2.6 

16% 

19% 

4.7 

26% 

30% 

4.9 

25% 

27% 

4.3 

22% 

Table  7.  R2  values  for  Image  Set  2 


Case 

R2 

R2  <  15  m 

2MSI 

0.8825 

0.8858 

3MSI 

0.9092 

0.8825 

2Pan 

0.8096 

0.8595 

3Pan 

0.8627 

0.8683 

46 


After  examining  the  multispectral  true  color  images  on  the  top  right  of  the  WKB 
output  figures  in  Appendices  A  and  B,  it  is  apparent  that  the  ocean  bottom  appears 
brighter  in  the  Image  Set  2  images.  This  is  most  likely  due  to  the  difference  in  collection 
angle  as  the  spacecraft  moved  over  the  area,  with  the  angle  becoming  more  favorable  for 
bottom  reflection  as  the  pass  evolved.  This  would  have  exacerbated  the  bottom  generated 
noise  issue  for  the  panchromatic  cases,  leading  to  the  further  degradation  in  performance. 

The  individual  depth  bin  values  tell  a  similar  story,  with  the  multispectral  cases 
doing  considerably  better  than  the  panchromatic  cases  through  a  depth  of  10  m,  where 
mean  depth  errors  were  nominally  1-2  m  better  and  relative  errors  were  10-40%  better  in 
the  multispectral  cases.  Above  10  m,  the  3Pan  case  performed  better  than  the  3MSI  case, 
with  mean  depth  errors  of  1.5-1. 7  m,  until  14  m  depth  where  3MSI  took  honors  again. 
Strangely,  the  multispectral  cases  did  not  degrade  as  rapidly  as  the  panchromatic  cases 
above  16  m,  unlike  in  set  one  where  they  degraded  roughly  equally. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

A  WKB  algorithm  was  investigated  to  determine  its  effectiveness  in  extracting 
bathymetry  from  WorldView-2  satellite  optical  imagery.  This  initial  analysis  indicates 
that  the  technique  described  herein  is  promising.  The  analysis  involved  comparing  WKB 
extracted  depths  to  SHOALS  LIDAR  survey  obtained  depths  at  Waimanalo  Beach  in 
Hawaii. 

The  algorithm  performed  well  in  the  midrange  depths  of  2-14  m,  its  best  cases 
managing  as  little  as  0.6  m  and  no  worse  than  2.3  m  of  mean  depth  error  in  any  depth  bin. 
Excluding  the  2-4  m  bin,  relative  errors  for  the  best  cases  ranged  from  a  best  of  10%  to  a 
worst  of  27%,  with  most  being  in  the  10-17%  range.  The  best  regression  line  slope 
achieved  was  1.01  with  an  associated  R~  of  91%. 

The  best  analysis  results  came  from  the  use  of  three  multispectral  images,  the 
worst  results  came  from  the  use  of  two  panchromatic  images.  The  panchromatic  cases 
would  certainly  have  performed  better  in  a  more  turbid  water  or  dark  bottom  area,  where 
the  bottom  would  not  have  contributed  to  a  poor  wave  SNR.  The  three-image  cases  did 
better  than  the  two-image  cases  in  all  but  a  couple  of  instances. 

There  may  have  been  some  error  introduced  into  the  comparisons  through  a 
conversion  of  the  original  imagery  from  latitude-longitude  geographic  coordinates  to 
UTM  coordinates,  which  may  have  differed  from  the  SHOALS  UTM  coordinates  by 
several  meters.  Understanding  that  the  bathymetry  hardly  changes  in  that  amount  of 
horizontal  distance,  and  that  the  WKB  algorithm  produced  depths  at  an  effective  spatial 
resolution  of  about  125  m,  this  error  should  be  insignificant.  Another,  perhaps  more 
significant,  source  of  error  is  the  elapsed  time  between  the  SHOALS  survey  and  the 
imagery  collection.  During  this  almost  11 -year  period,  the  bathymetry  may  have 
changed,  perhaps  significantly  in  some  portions. 

In  spite  of  these  potential  sources  of  error  and  given  the  fairly  complex 
environment  in  which  it  was  tested,  the  performance  of  the  algorithm  was  notable. 
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Unlike  many  previous  WKB  studies,  which  were  conducted  in  areas  with  simpler,  nearly- 
constant  sloping  bottoms,  often  assumed  cross-shore  currents,  and  avoided  wave 
refraction,  Waimanalo  Beach  has  fairly  complex  bottom  topography,  intricate  surface 
currents,  and  noticeable  wave  refraction. 

WKB  is  limited  in  its  ability  to  do  depth  extraction  in  the  surf  zone  (less  than 
about  2  m)  and  beyond  some  deep  limit  determined  by  the  surface  wave  spectrum.  In 
this  case,  the  lack  of  long  wave  swells  limited  WKB’s  ability  to  determine  depth 
satisfactorily  to  approximately  15  m.  A  scene  with  more  long  wave  swells  in  it  probably 
would  have  produced  more  accurate  depths  out  to  20  m  or  more. 

B.  RECOMMENDATIONS  FOR  FUTURE  WORK 

The  WKB  algorithms  and  codes  are  still  evolving.  Continuing  work  will 
determine  whether  these  changes  improve  the  algorithm’s  results.  Among  the  newly 
implemented  updates  is  the  fusion  of  panchromatic  and  multispectral  images  into  the 
same  depth  extraction  routine,  whereas  they  were  calculated  separately  before.  This 
allows  the  algorithm  to  take  advantage  of  the  strengths  of  each  image  set,  using  the 
higher  resolution  panchromatic  images  for  the  current  determination  and  the  bands  of  the 
multispectral  images  that  create  the  best  SNR  to  determine  the  depth. 

The  code  based  on  the  tiling  approach,  used  here,  was  recently  replaced  with  code 
that  makes  the  algorithm  more  computationally  efficient  with  the  non-tiling  method.  One 
way  to  really  capitalize  on  this  is  to  introduce  parallel  processing  into  the  process.  This 
would  significantly  speed  up  the  computing  time  for  larger  areas  and  more  image  pairs. 

Future  improvements  include:  a  wave  contrast  metric  to  help  identify  when 
contrast  has  played  a  role  in  the  results;  better  logic  to  improve  performance  with 
multiple  image  pairs;  a  way  to  flag  extremes  and  replace  them  with  appropriate  values  to 
eliminate  unreasonable  depth  conclusions  in  the  surf  zone  or  at  the  deep  limits;  and 
accuracy  and  spatial  resolution  improvements. 

In  addition  to  testing  algorithm  improvements,  future  work  could  investigate 
different  combinations  of  multispectral  bands  to  see  which  produces  the  best  results  by 
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better  isolating  the  waves  from  the  background.  Perhaps  principle  component  transforms 
of  the  images  can  be  applied  to  increase  wave  visibility  and  SNR. 

Improvements  in  data  collection  could  also  be  explored.  For  example,  reducing 
the  time  interval  between  subsequent  image  pairs  would  reduce  wave  propagation  and  the 
change  in  view  angle  from  one  image  to  the  next.  This  will  help  with  accuracy  by 
reducing  the  opportunity  for  the  waves  to  be  modified  by  sources  other  than  depth 
(making  the  dispersion  relation  a  more  exact  predictor  of  the  wave  motion)  and  by 
allowing  the  multiple-image-pair  averaging  technique  to  have  more  success.  It  would 
also  alleviate  the  ambiguity  issue  associated  with  sampling  at  the  same  time  interval  as 
the  period  of  a  bulk  of  the  waves. 

There  are  other  ways  to  realize  the  last  two  gains  if  reducing  the  image  collection 
time  interval  becomes  too  difficult.  By  incorporating  a  view-angle  transfer  function  to 
account  for  the  change  in  angle,  the  averaging  technique  could  still  produce  improvement 
when  averaging  the  image  pairs.  The  ambiguity  problem  can  probably  be  overcome  with 
a  better  algorithm  to  more  effectively  deal  with  it. 

Finally,  better  choices  in  the  imagery  product  could  improve  future  efforts. 
Receiving  the  imagery  in  the  UTM  format  that  WKB  requires  would  eliminate  the  need 
to  convert  it,  avoiding  any  error  introduced  in  the  process.  Obtaining  accurate  timing 
details  from  the  vendor  to  establish  more  precisely  when  each  image  line  was  scanned 
will  remove  the  need  to  estimate  during  the  pre-processing.  While  the  estimation  is  good, 
knowing  the  timing  exactly  will  help  improve  WKB  accuracy. 
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Figure  25.  2MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  26.  Comparison  maps  of  SFIOALS  (left)  and  2MSI  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


Figure  27.  2MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 

just  depths  less  than  15  m  (bottom) 
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WKB  Depth  Error  vs  Depth 


Depth  (m) 

Figure  28.  2MSI  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  29.  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  30.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  3 1 .  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  32.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  33.  3MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  34.  Comparison  maps  of  SHOALS  (left)  and  3MSI  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


3MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 
just  depths  less  than  15  m  (bottom) 
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Figure  35. 


WKB  Depth  Error  vs  Depth 


Depth  (m) 

Figure  36.  3MSI  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  37.  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  38.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  39.  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  40.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  4 1 .  2Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  42.  Comparison  maps  of  SHOALS  (left)  and  2Pan  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


Figure  43.  2Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 

just  depths  less  than  15  m  (bottom) 
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WKB  Depth  Error  vs  Depth 


Depth  (m) 

Figure  44.  2Pan  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  45.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  46.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  47.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  48.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  49.  3Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  50.  Comparison  maps  of  SHOALS  (left)  and  3Pan  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


Figure  5 1 .  3Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 

just  depths  less  than  15  m  (bottom) 
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WKB  Depth  Error  vs  Depth 


Depth  (m) 

Figure  52.  3Pan  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  53.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  54.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  55.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  56.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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APPENDIX  B.  IMAGE  SET  2 
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Figure  57.  2MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  58.  Comparison  maps  of  SFIOALS  (left)  and  2MSI  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


2MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 
just  depths  less  than  15  m  (bottom) 
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Figure  59. 


WKB  Depth  Error  vs  Depth 


Depth  (m) 

Figure  60.  2MSI  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  6 1 .  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  62.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  63.  Transect  at  specified  location  showing  2MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 


WKB-SHOALS  Depth  Comparison  at  2363  UTM  Northing  (km) 


WKB  Relative  Depth  Error  at  2363  UTM  Northing  (km) 


UTM  Easting  (km) 


Figure  64.  Transect  at  specified  location  showing  2MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  65.  3MSI  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  66.  Comparison  maps  of  SFIOALS  (left)  and  3MSI  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


3MSI  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 
just  depths  less  than  15  m  (bottom) 
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Figure  67. 


WKB  Depth  Error  vs  Depth 
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Figure  68.  3MSI  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  69.  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  70.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  7 1 .  Transect  at  specified  location  showing  3MSI  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  72.  Transect  at  specified  location  showing  3MSI  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  73.  2Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  74.  Comparison  maps  of  SHOALS  (left)  and  2Pan  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


Figure  75.  2Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 

just  depths  less  than  15  m  (bottom) 
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Figure  76.  2Pan  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  77.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  78.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  79.  Transect  at  specified  location  showing  2Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  80.  Transect  at  specified  location  showing  2Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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Figure  8 1 .  3Pan  WKB  output  showing  wave  direction  (top  left),  a  true  color 
image  (top  right),  extracted  bathymetry  (bottom  left),  and  extracted  ocean 

currents  (bottom  right) 
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Figure  82.  Comparison  maps  of  SHOALS  (left)  and  3Pan  WKB  (right) 


WKB-SHOALS  Comparison 


WKB-SHOALS  Comparison  with  a  cutoff  of -15  m 


Figure  83.  3Pan  scatter  plot  showing  thinned  data  for  all  depths  (top)  and  for 

just  depths  less  than  15  m  (bottom) 
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Mean  Depth  Error  (m) 
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WKB  Depth  Error  vs  Depth 


Figure  84.  3Pan  bar  graph  showing  the  mean  depth  error  for  several  depth 

bins 
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Figure  85.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 
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Figure  86.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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WKB-SHOALS  Depth  Comparison  at  2362  UTM  Northing  (km) 


Figure  87.  Transect  at  specified  location  showing  3Pan  WKB  (see  below  for 

more  explanation  of  the  panels) 


WKB-SHOALS  Depth  Comparison  at  2363  UTM  Northing  (km) 
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Figure  88.  Transect  at  specified  location  showing  3Pan  WKB  depth  plotted 
with  SHOALS  depth  (top);  the  difference  between  them,  or  depth  error 
(mid);  and  depth  error  as  a  percentage  of  depth,  or  relative  depth  error  (hot) 
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